home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / ssymm.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  9.7 KB  |  298 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE SSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
  5.      $                   BETA, C, LDC )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        SIDE, UPLO
  8.       INTEGER            M, N, LDA, LDB, LDC
  9.       REAL               ALPHA, BETA
  10. *     .. Array Arguments ..
  11.       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  SSYMM  performs one of the matrix-matrix operations
  18. *
  19. *     C := alpha*A*B + beta*C,
  20. *
  21. *  or
  22. *
  23. *     C := alpha*B*A + beta*C,
  24. *
  25. *  where alpha and beta are scalars,  A is a symmetric matrix and  B and
  26. *  C are  m by n matrices.
  27. *
  28. *  Parameters
  29. *  ==========
  30. *
  31. *  SIDE   - CHARACTER*1.
  32. *           On entry,  SIDE  specifies whether  the  symmetric matrix  A
  33. *           appears on the  left or right  in the  operation as follows:
  34. *
  35. *              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
  36. *
  37. *              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
  38. *
  39. *           Unchanged on exit.
  40. *
  41. *  UPLO   - CHARACTER*1.
  42. *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  43. *           triangular  part  of  the  symmetric  matrix   A  is  to  be
  44. *           referenced as follows:
  45. *
  46. *              UPLO = 'U' or 'u'   Only the upper triangular part of the
  47. *                                  symmetric matrix is to be referenced.
  48. *
  49. *              UPLO = 'L' or 'l'   Only the lower triangular part of the
  50. *                                  symmetric matrix is to be referenced.
  51. *
  52. *           Unchanged on exit.
  53. *
  54. *  M      - INTEGER.
  55. *           On entry,  M  specifies the number of rows of the matrix  C.
  56. *           M  must be at least zero.
  57. *           Unchanged on exit.
  58. *
  59. *  N      - INTEGER.
  60. *           On entry, N specifies the number of columns of the matrix C.
  61. *           N  must be at least zero.
  62. *           Unchanged on exit.
  63. *
  64. *  ALPHA  - REAL            .
  65. *           On entry, ALPHA specifies the scalar alpha.
  66. *           Unchanged on exit.
  67. *
  68. *  A      - REAL             array of DIMENSION ( LDA, ka ), where ka is
  69. *           m  when  SIDE = 'L' or 'l'  and is  n otherwise.
  70. *           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
  71. *           the array  A  must contain the  symmetric matrix,  such that
  72. *           when  UPLO = 'U' or 'u', the leading m by m upper triangular
  73. *           part of the array  A  must contain the upper triangular part
  74. *           of the  symmetric matrix and the  strictly  lower triangular
  75. *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  76. *           the leading  m by m  lower triangular part  of the  array  A
  77. *           must  contain  the  lower triangular part  of the  symmetric
  78. *           matrix and the  strictly upper triangular part of  A  is not
  79. *           referenced.
  80. *           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
  81. *           the array  A  must contain the  symmetric matrix,  such that
  82. *           when  UPLO = 'U' or 'u', the leading n by n upper triangular
  83. *           part of the array  A  must contain the upper triangular part
  84. *           of the  symmetric matrix and the  strictly  lower triangular
  85. *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  86. *           the leading  n by n  lower triangular part  of the  array  A
  87. *           must  contain  the  lower triangular part  of the  symmetric
  88. *           matrix and the  strictly upper triangular part of  A  is not
  89. *           referenced.
  90. *           Unchanged on exit.
  91. *
  92. *  LDA    - INTEGER.
  93. *           On entry, LDA specifies the first dimension of A as declared
  94. *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  95. *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
  96. *           least  max( 1, n ).
  97. *           Unchanged on exit.
  98. *
  99. *  B      - REAL             array of DIMENSION ( LDB, n ).
  100. *           Before entry, the leading  m by n part of the array  B  must
  101. *           contain the matrix B.
  102. *           Unchanged on exit.
  103. *
  104. *  LDB    - INTEGER.
  105. *           On entry, LDB specifies the first dimension of B as declared
  106. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  107. *           max( 1, m ).
  108. *           Unchanged on exit.
  109. *
  110. *  BETA   - REAL            .
  111. *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
  112. *           supplied as zero then C need not be set on input.
  113. *           Unchanged on exit.
  114. *
  115. *  C      - REAL             array of DIMENSION ( LDC, n ).
  116. *           Before entry, the leading  m by n  part of the array  C must
  117. *           contain the matrix  C,  except when  beta  is zero, in which
  118. *           case C need not be set on entry.
  119. *           On exit, the array  C  is overwritten by the  m by n updated
  120. *           matrix.
  121. *
  122. *  LDC    - INTEGER.
  123. *           On entry, LDC specifies the first dimension of C as declared
  124. *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  125. *           max( 1, m ).
  126. *           Unchanged on exit.
  127. *
  128. *
  129. *  Level 3 Blas routine.
  130. *
  131. *  -- Written on 8-February-1989.
  132. *     Jack Dongarra, Argonne National Laboratory.
  133. *     Iain Duff, AERE Harwell.
  134. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  135. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  136. *
  137. *
  138. *     .. External Functions ..
  139.       LOGICAL            LSAME
  140.       EXTERNAL           LSAME
  141. *     .. External Subroutines ..
  142.       EXTERNAL           XERBLA
  143. *     .. Intrinsic Functions ..
  144.       INTRINSIC          MAX
  145. *     .. Local Scalars ..
  146.       LOGICAL            UPPER
  147.       INTEGER            I, INFO, J, K, NROWA
  148.       REAL               TEMP1, TEMP2
  149. *     .. Parameters ..
  150.       REAL               ONE         , ZERO
  151.       PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  152. *     ..
  153. *     .. Executable Statements ..
  154. *
  155. *     Set NROWA as the number of rows of A.
  156. *
  157.       IF( LSAME( SIDE, 'L' ) )THEN
  158.          NROWA = M
  159.       ELSE
  160.          NROWA = N
  161.       END IF
  162.       UPPER = LSAME( UPLO, 'U' )
  163. *
  164. *     Test the input parameters.
  165. *
  166.       INFO = 0
  167.       IF(      ( .NOT.LSAME( SIDE, 'L' ) ).AND.
  168.      $         ( .NOT.LSAME( SIDE, 'R' ) )      )THEN
  169.          INFO = 1
  170.       ELSE IF( ( .NOT.UPPER              ).AND.
  171.      $         ( .NOT.LSAME( UPLO, 'L' ) )      )THEN
  172.          INFO = 2
  173.       ELSE IF( M  .LT.0               )THEN
  174.          INFO = 3
  175.       ELSE IF( N  .LT.0               )THEN
  176.          INFO = 4
  177.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  178.          INFO = 7
  179.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  180.          INFO = 9
  181.       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
  182.          INFO = 12
  183.       END IF
  184.       IF( INFO.NE.0 )THEN
  185.          CALL XERBLA( 'SSYMM ', INFO )
  186.          RETURN
  187.       END IF
  188. *
  189. *     Quick return if possible.
  190. *
  191.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  192.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  193.      $   RETURN
  194. *
  195. *     And when  alpha.eq.zero.
  196. *
  197.       IF( ALPHA.EQ.ZERO )THEN
  198.          IF( BETA.EQ.ZERO )THEN
  199.             DO 20, J = 1, N
  200.                DO 10, I = 1, M
  201.                   C( I, J ) = ZERO
  202.    10          CONTINUE
  203.    20       CONTINUE
  204.          ELSE
  205.             DO 40, J = 1, N
  206.                DO 30, I = 1, M
  207.                   C( I, J ) = BETA*C( I, J )
  208.    30          CONTINUE
  209.    40       CONTINUE
  210.          END IF
  211.          RETURN
  212.       END IF
  213. *
  214. *     Start the operations.
  215. *
  216.       IF( LSAME( SIDE, 'L' ) )THEN
  217. *
  218. *        Form  C := alpha*A*B + beta*C.
  219. *
  220.          IF( UPPER )THEN
  221.             DO 70, J = 1, N
  222.                DO 60, I = 1, M
  223.                   TEMP1 = ALPHA*B( I, J )
  224.                   TEMP2 = ZERO
  225.                   DO 50, K = 1, I - 1
  226.                      C( K, J ) = C( K, J ) + TEMP1    *A( K, I )
  227.                      TEMP2     = TEMP2     + B( K, J )*A( K, I )
  228.    50             CONTINUE
  229.                   IF( BETA.EQ.ZERO )THEN
  230.                      C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
  231.                   ELSE
  232.                      C( I, J ) = BETA *C( I, J ) +
  233.      $                           TEMP1*A( I, I ) + ALPHA*TEMP2
  234.                   END IF
  235.    60          CONTINUE
  236.    70       CONTINUE
  237.          ELSE
  238.             DO 100, J = 1, N
  239.                DO 90, I = M, 1, -1
  240.                   TEMP1 = ALPHA*B( I, J )
  241.                   TEMP2 = ZERO
  242.                   DO 80, K = I + 1, M
  243.                      C( K, J ) = C( K, J ) + TEMP1    *A( K, I )
  244.                      TEMP2     = TEMP2     + B( K, J )*A( K, I )
  245.    80             CONTINUE
  246.                   IF( BETA.EQ.ZERO )THEN
  247.                      C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2
  248.                   ELSE
  249.                      C( I, J ) = BETA *C( I, J ) +
  250.      $                           TEMP1*A( I, I ) + ALPHA*TEMP2
  251.                   END IF
  252.    90          CONTINUE
  253.   100       CONTINUE
  254.          END IF
  255.       ELSE
  256. *
  257. *        Form  C := alpha*B*A + beta*C.
  258. *
  259.          DO 170, J = 1, N
  260.             TEMP1 = ALPHA*A( J, J )
  261.             IF( BETA.EQ.ZERO )THEN
  262.                DO 110, I = 1, M
  263.                   C( I, J ) = TEMP1*B( I, J )
  264.   110          CONTINUE
  265.             ELSE
  266.                DO 120, I = 1, M
  267.                   C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
  268.   120          CONTINUE
  269.             END IF
  270.             DO 140, K = 1, J - 1
  271.                IF( UPPER )THEN
  272.                   TEMP1 = ALPHA*A( K, J )
  273.                ELSE
  274.                   TEMP1 = ALPHA*A( J, K )
  275.                END IF
  276.                DO 130, I = 1, M
  277.                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
  278.   130          CONTINUE
  279.   140       CONTINUE
  280.             DO 160, K = J + 1, N
  281.                IF( UPPER )THEN
  282.                   TEMP1 = ALPHA*A( J, K )
  283.                ELSE
  284.                   TEMP1 = ALPHA*A( K, J )
  285.                END IF
  286.                DO 150, I = 1, M
  287.                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
  288.   150          CONTINUE
  289.   160       CONTINUE
  290.   170    CONTINUE
  291.       END IF
  292. *
  293.       RETURN
  294. *
  295. *     End of SSYMM .
  296. *
  297.       END
  298.